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Abstract. We analyze the flow of a liquid foam between two plates separated by a gap of the order of 
the bubble size (2D foam). We concentrate on the salient features of the flow that are induced by the 
presence, in an otherwise monodisperse foam, of a single large bubble whose size is one order of magnitude 
larger than the average size. We describe a model suited for numerical simulations of flows of 2D foams 
made up of a large number of bubbles. The numerical results are successfully compared to analytical 
predictions based on scaling arguments and on continuum medium approximations. When the foam is 
pushed inside the cell at a controlled rate, two basically different regimes occur: a plug flow is observed 
at low flux whereas, above a threshold, the large bubble migrates faster than the mean flow. The detailed 
characterization of the relative velocity of the large bubble is the essential aim of the present paper. The 
relative velocity values, predicted both from numerical and from analytical calculations that are discussed 
here in great detail, are found to be in fair agreement with experimental results from pQ. 
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1 Introduction each individual cell of a 2D foam can be followed during 

its motion and its deformation. The relation between the 

* , n . , , -ln-jr local structure and the induced macroscopic visco-elastic 

Among complex fluids and structured materials, fluid loams, 

, . ., • i r j- j • i flow can thus be directly observed. These reasons explain 

which are primarily materials oi outstanding industrial 

, "i , , why theoretical studies of 2D foam flows in viscous regimes 

importance, deserve too a particular attention as model ■' ° 

, T j , , ii i4-* i u are of growing importance 

systems. Indeed, they are governed by relatively well un- ° ° 1 

derstood and simple local equilibrium properties. Further, 
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satisfactory description of the underlying elementary in- 
stability which involves a single defect. Preliminary results 
on this instability were shortly discussed in a previous arti- 
cle jS] , and experimental results are thoroughly described 
in £Q. In the following, we shall report in detail on the- 
oretical and numerical investigations of the flow of a 2D 
monodisperse foam in which a single large bubble is em- 
Fig. 1. Schematic drawing of the set-up from which we col- bedded, leading to what is called hereafter the large bubble 
lected the experimental results 0f referred to in the present instability 
theoretical paper. 
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By a 2D foam, we refer here to foams confined between 
two horizontal plates and made up of a single bubble layer. 
The richness of the flow properties is largely related to the 
disorder of the foam structure ^U]> but we shall show in 
this paper that even few defects in an otherwise basically 
ordered structure may already induce strong modifications 
of the flow in a viscous regime. 

If a constant foam flux is imposed in a Hele-Shaw cell 
(see Fig^), a uniform pressure gradient is created between 
upstream and downstream to compensate the viscous drag 
on the plates. Lateral boundary effects are negligible and 
the reference flow at low velocity is thus a simple plug 
flow. At higher velocity, the plug flow becomes unstable 
for polydisperse foams. The velocity field, then dominated 
by the bubble size distribution and by its spatial correla- 
tions, leads to a very complex and fluctuating flow. Larger 
bubbles move faster than the smaller ones and the foam 
structure is fully reorganized with strong attractive corre- 
lations between bubbles of similar sizes. An understanding 
of this surprising foam behavior is out of reach without a 



The viscous dissipation is localized on the Plateau bor- 
ders touching the plates, ie the contact lines between the 
vertical foam films and the films wetting the plates. In a 
monodisperse foam which flows uniformly, viscous forces, 
averaged at the scale of few small bubbles, induce a uni- 
form stress field oriented upstream. For a large bubble 
(denoted by LB in the following) the density of Plateau 
borders is smaller, and the averaged viscous stress is lower. 
The LB is, in a way, analogous to a low density drop em- 
bedded in a fluid submitted to a uniform gravity field. In 
both cases, the uniform stress field is balanced by a lin- 
ear pressure field, and a resulting Archimedes-like force 
is exerted on the bubble (or respectively on the drop) 
in the direction opposite to that of the force field. This 
driving force, responsible for the LB migration, competes 
with the elastic reaction of the bubble network : a stable 
foam structure is obtained at low velocity. By contrast, 
the large bubble starts to migrate through the small bub- 
ble network, faster than the mean flow, when the foam 
velocity is large enough to induce a driving effect on the 
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large bubble greater larger than the plastic threshold of which the foam is described by the positions of its vertices, 



the foam (see Fig. 01. 




Fig. 2. Two 'shots' of the LB migration through the foam, in 
the LB frame (numerical simulations). The foam crystalline 
orientation is the most favorable and the migration thus occurs 
in the same direction as the mean flow, to the right. 

Our main theoretical results consist of two scaling laws 
derived from general dimensional arguments for the ve- 
locity threshold (eq. L'MI) and for the large bubble velocity 
beyond that threshold (eq. . They are compared with 
the results of numerical simulations based on the vertex 
model ||)| and with those of experiments presented in pQ. 

The paper is organized as follows. In Section |21 we 
discuss the main ingredients (bubble shape, viscous force, 
etc.) that play a role in the process by which the insta- 
bility takes place. Section [3] presents a general frame, in 



that may be used for numerical simulations of any 2D 
foam flow. Section 21 is devoted to the qualitative descrip- 
tion of the performed simulations and of the flow behavior 
so obtained. An analytical model, based on a continuous 
description, is introduced in Section [5] and allows quan- 
titative comparison with the numerical results described 
in Section |SJ The pressure field in the whole foam is dis- 
cussed in Section [7| Finally, the results and perspectives 
are summarized in Section [S] 

2 Specific properties of the 2D foam. 
2.1 Equilibrium film shape. 

The two-dimensional foam referred to here is a foam con- 
fined between to plates separated by a distance h, much 
smaller than the typical bubble size. For small values of h, 
the bubbles are organized in a single layer and every film 
touches the upper and the lower plates The Plateau 
borders touching the upper (or lower) plate draw a net- 
work r of connected curved lines, as shown on FigFJJ For 
2D foams in equilibrium, the films separating two bubbles 
are part of vertical cylinders and the only non vanishing 
curvature is thus the constant curvature in a plane paral- 
lel to the plates, denoted by C\ in the following. Thus P 
contains the whole information on the foam structure. 

This property entails many important consequences 
specific to 2D foams in equilibrium 12^ and yields a way 
of measuring the bubble pressure by simple image analy- 
sis using the Laplace relation AP — jci , where AP is the 
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pressure difference between two adjacent bubbles and 7 play here the role of the liquid reservoir and the relative 
the surface tension. In a flowing foam, the viscous forces, motion between the liquid and the solid plate is due to the 
which are localized close to the plates, modify the con- motion of the bubble walls. According to the Bretherton 
tact angle 9 as depicted on Fig|31and0| Both curvatures theory [^Oj, the highest velocity gradients occur in a small 
may be then of the same order of magnitude as evaluated domain connecting the Plateau borders and the wetting 
in section 12.31 and the pressure can no more be read on films, assumed to be at rest on each plate. This region is 
a foam image, as the mean curvature is not measurable of characteristic thickness 21) 
easily anymore. 

6{v) ~ Rpii^vh) 2 ' 3 -510- 7 ™ (1) 

2.2 Viscous forces. with r/ u , = 10~ 3 Pa s the liquid phase viscosity, 7 ~ 30 10~ 3 N/m 

the surface tension, v ~ 10~ 2 m/s the macroscopic velocity 
of the Plateau border and Rpi ~ 10~ 4 m the Plateau bor- 
der size. The relevant Reynolds number for this problem 
is thus Re — p5v/r) w ~ 5 10~ 2 , whose magnitude allows 
us to neglect the non-linear convective term in the Navier- 
Stokes equation. The viscous force /° (per unit length of 
Plateau border) exerted on a Plateau border sliding at 
a velocity v — vu v on the plate varies nevertheless non 
linearly with v, namely as a power law with an exponent 
smaller than unity. It will be written /° = —rj{v)v, rj(v) 
being an effective viscosity (with the dimension of a dy- 
namic vioscosity). As the films swell when the foam ve- 
locity increases (see eq. 0, velocity variations occur on 
larger length scales and the effective viscosity decreases. 
We assume here that r](v) ~ v^ 1 / 3 (Landau Levich expo- 
nent), although power laws with exponent values ranging 
between —1/3 and —1/2 are reported in the literature |17| . 
This force per unit length of Plateau border is finally of 
the order of 10~ 3 N/m and is expressed as ^5] : 



The effective viscosity of a liquid foam pushed between 
two plates, ie the driving force to velocity ratio, is typi- 
cally one or two orders of magnitude larger than that of 
the pure liquid phase, despite the very large amount of 
gas (> 95%) trapped in the material [2lll3U14U15| . This 
surprising behavior is due to the existence of very small 
length scales at which the liquid phase is confined in the 
structure. In the 2D geometry, the dominating process is 
the force between the Plateau borders and the plates. The 
viscous dissipation may have various origins, including the 
bulk viscosity of the liquid soap water, the viscosity of the 
surfactant monolayer and the diffusion resistance, leading 
to different relations between pressure drop and velocity 
[161117) . Considering a fluid surfactant monolayer, we as- 
sume here that the dominating contribution comes from 
the bulk viscosity. The case of a rigid monolayer will nev- 
ertheless be considered too in section I5~2l 

The local flow geometry close to the plates is very simi- 
lar to the well-known Landau-Levitch situation, in which a 
plate is pulled out of a liquid |18II19| . The Plateau borders /° = —rj(v)v = — XrjwvCa^^^y ■ n\u v ; . (2) 
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with Ca — r) w v/d the capillary number of the order of 
3 10~ 4 . The numerical prefactor A ~ 10 and the geomet- 
rical factor were measured in similar flow conditions, the 
latter involving u v the unit vector of the film velocity di- 
rection and n the unit vector normal to the film (with 
an arbitrary orientation) 15 |. The force direction is not 
clearly evidenced and seems difficult to determine directly. 
It is here assumed to be in the direction opposite to that 
of the film velocity. 

It will be useful to express instead the viscous force 
per unit surface of the 2D foam averaged on the bubble 
scale (mean viscous stress on the plates). With a typical 
bubble size of d ~ 10~ 2 m we obtain 



02 oi 
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Ca~ 1/3 u v ~ W^N/m 2 



(3) 



For comparison the viscous stress obtained for pure water 
with the same mean velocity is of the order of rj w v/h ~ 
lQ- 2 N/m 2 . 

2.3 Out of equilibrium film shape 




Fig. 3. A schematic side view of the Plateau border during 
motion. 




Fig. 4. The film curvatures. The Plateau border curvature is 
locally ci = 1/Ri- In the plane normal to the plates and to the 
Plateau border the film section can be approximated by an arc 
of circle of radius R2 ■ 

During flow, as inertia and dissipation in the bubble 
walls are negligible, bubbles still obey the Laplace equa- 
tion AP = 7c where c is the sum of the two main curva- 
tures which remains constant on the whole bubble wall. 
The boundary conditions on the plates are modified by 
the viscous force, leading to another film shape. The force 
balance on an infinitesimal volume around a Plateau bor- 
der (see Fig. |3J) allows to determine the orientation 9 of 
the fluid film at the contact with the wall, 



7sin0 = /° 



(4) 



The film shape is thus a surface of constant mean curva- 
ture with a contact angle at the plates depending on the 
local velocity. The two principal curvatures vary on the 
surface but may be approximated by c p \ ~ c\ = 1/ R\ and 
c P 2 ~ C2 = I/-R2, with ci the curvature in the symmetry 
plane between the two plates and C2 = 2sm9/h ~ iOm^ 1 
(see Fig. 0J. Both curvatures may thus be of the same 
order of magnitude for centimeter-sized bubbles. As the 
mean curvature remains constant on a given film, if ci 
varies along the film, then c\ varies too, and the Plateau 
border is no more an arc of circle. 
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3 Numerical simulation of a 2D flow 
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Simulations or analytical calculations make a compromise 
between the complexity of the dissipative function they 
use [SHE] and the number of bubbles they can consider 
r4||6||2||8]|9] . The numerical model presented in this section 
is well adapted for a broad family of flow configurations 
in 2D. It is based on a realistic dynamical behavior, even Fig. 5. Force balance on the system Si around the vertex i 




at large velocities, and on a simplified foam structure de- 
scription which allows to use a large number of bubbles at 
reasonable computational cost. 

3.1 Numerical variables and equation of motion 

As discussed above, the 3D film shapes are unknown. A 
local description of the 2D Plateau borders network r 
has been used by Kern et al. in numerical simulations 
performed with a small number of bubbles in which only 



is then defined as being composed of the vertex i with its 
three outgoing edges reduced by half, as depicted on Fig. [5] 
The equation of motion is provided by the force balance on 
each system Si. The number of vertices remains constant 
as we forbid the occurrence of film breakage and of foam 
coarsening. 

The tension force is, with r,j = — Tj and = 



7^E? (5) 



36Ji T%3 



c\ is taken into account while C2 [H| is neglected. This ap- We define to be the normal to the edge (ij), oriented 

proach requires to discretize every Plateau border line and arbitrarily, say from a bubble A towards a bubble B, and 

becomes very demanding for large scale simulations. For SPij(— Pb — Pa) to be the pressure jump on this edge, 

numerical efficiency, we choose to simplify the structure The resulting pressure force is then : 
further by assuming vertical films to be planar, as already „ 
done in The Plateau borders touching the plates are jeJi 

straight lines (called edges), entirely determined by their Finally the viscous force determination imposes some 

end points (the vertices) . Our numerical results have been additional approximations. The whole system S 1 , is as- 

quantitatively reproduced by Sanyal et al. with a method sumed to move with the same velocity as the vertex (i) 

based on a Potts model [2T]. : v » = v i u i- Thus we get 
The variables used are the vertex locations r,. The 



k Yl \\ ri i x u *k 2/3u * ( 7 ) 

points located at and Tj with j 6 J^. The system Si with the prefactor k obtained from eq. k = Arfw ; 7 1 / 3 . 



2 

three edges which meet at a given vertex i have their end i&Ji 
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The equation of motion, F v ,i + F t ^ + F p .i = im- the same for every bubble in the foam, in contrast to the 
poses the orientation of Vi along the unit vector Fi/Fi, local reference P eq that will be introduced in the following 
with Fi = Ffj + Fp^. We get thus the following explicit paragraphs. In our simulations, the reference pressure P° 
expression for the velocity, as a function of the positions is strongly underestimated to enhance the numerical sta- 
rj, and of Fi which can itself be expressed as a function bility. Nevertheless, even if P° is smaller than its experi- 
of the positions and the pressures fea l5l6fl : mental value, it remains much larger than all the pressure 

variations induced by the flow and the system is in the 
v — F F 2 I k | — r • • x F | J (8) incompressible limit (the relative variations of the bubble 

areas remain smaller than few %). The pressure field given 

As we will see now, the pressures are also related to the 

by the product P°SA/A is thus independent on P°, even 

positions, so that the velocity is a function of the only 

if both factors obviously depend on it. The choice of the 

variables Tj. An iterative process is then possible as, for 

compressibility value (or equivalently of the pressure ref- 

each time step, the velocities and consequently the dis- 

erence) has therefore no real influence on the dynamical 

placements of the vertices can be computed, leading to an 

behavior, as checked numerically. 

actualisation of their positions. 

This formalism presents the big advantage to provide 
an explicit expression of the pressure as a function of the 

3.2 Explicit expression of the pressure 

bubble geometry, in contrast to the classical Lagrange 
On the time scale of a few seconds considered in this pa- multiplier approach in which a strict incompressibility is 
per, gas diffusion is negligible but pressure equilibrium is assumed. 

reached. The pressure is thus uniform in each bubble k and The rererence area is chosen for each bubble as a fixed 
can be determined, under the assumption of an ideal gas, parameter . The actual area is computed after each dis- 
as a function of the area A k of the bubble in contact with placement as a f uncti on of the positions of the vertices, 
the plate, knowing the constant quantity of gas n k trapped assuming a po lyhedral shape for the bubble. Finally, the 
inside and its isothermal compressibility X t = 1/Pq ■ pressure value in the bubble is deduced from eq. E 

p k = — ^ fc '° r~j p° _ pO (^fc ~ ^*ig) ^g-j Note that the Plateau rules which constrain the values 

Ak Ak,o 

of the angles between the tangents of connected films at 

where 

vertices, which are for instance equal to 120° in a foam 

A k o = n k R PG T/hP° (10) 

at equilibrium, are not imposed in this modeling because 
is a reference value, obtained if the bubble pressure is the films are represented by straight lines. In this case, the 
equal to the global reference pressure P°. This pressure is pressure difference between adjacent bubbles is no more 
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related to a film curvature but to a deviation of the angles 
from the reference value of 120°. 

3.3 Building of an initial foam structure 

Periodic boundary conditions are imposed to the foam 
structure on the sides of the rectangular simulation box, 
the direction of the flow being parallel to the longer side. 
First of all, an artificial network of a few thousands of 
bubbles is built, with the positions and the connectivities 
of all vertices. Most results were obtained with an initial 
structure built from a perfect hexagonal network, but few 
simulations were performed with a disordered foam ob- 
tained from a Voronoi tessellation. Then a reference area 
is attributed to each bubble to reach the sought-after foam 
polydispersity and to create possibly very large bubbles 
in the foam. The system Sq is then relaxed to an equilib- 
rium foam structure, using the iterative process detailed 
in paragraph 13.41 until the vertex motions become of the 
order of the numerical noise. After this first relaxation 
step, the areas are very close to the reference areas Akfi. 

3.4 Temporal evolution 

The foam evolution is based on eq|SJ It is determined it- 
eratively in time from an initial structure Sq , that may be 
cither an equilibrium shape or an arbitrarily chosen out-of- 
equilibrium shape. Both transient or stationary processes 
may then be studied either in the presence of external 
perturbations or during relaxation to equilibrium. If we 
want to impose a flow in the direction x, we choose a set 
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of vertices connected by a wiggly line C of edges whose 
average orientation is along y (Fig.^J). For that purpose, 
we select first a vertex with a first neighbor on the other 
side of the line y — 0. Then, from this point n, we select 
the neighbor of smallest |x| verifying y > j/j, and we it- 
erate the process until we recover the initial point. The 
motion of the vertices which belong to that set C does not 
obey eq. |H1 as do the other vertices, but they are forced 
to move at the desired velocity at each time step. This 
induces an artificial pressure discontinuity across the line 
C, connecting the upstream and the downstream regions 
in the periodic box. The information propagates through 
the whole foam thanks to the pressure field. 

The main algorithm is the following : During the time 
step dt, (a) the vertices belonging to £ are displaced from 
Vodi(m case of flow), (b) the bubble areas and pressures 
are updated, (c) Velocities, and thus displacements, are 
computed everywhere else with eq|Hl (d) If the vertex dis- 
placements result in the formation of too small edges, Tl 
events are performed (see below), otherwise, (d') the dis- 
placements vdt are performed. 

If an edge becomes smaller than a fixed value e, a topo- 
logical transformation Tl is performed, as explained on 
Fig- El F° r stability reasons, only one Tl is allowed to oc- 
cur at each time step in the whole foam. In the rare case 
where two edges become simultaneously too small, we sim- 
ply choose to perform only one of the two transformations, 
the other being done at the next time step. This would 
have to be improved to deal with different situations, for 
instance the case of foam shearing. The edge is switched 
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Fig. 6. The encircled vertices belong to the line C. Their ve- 
locities are fixed to a given value to impose the flow rate. The 
two axes are sides of the simulation box, gray vertices result 
from the use of periodic boundary conditions. 




Fig. 7. The structure before (dashed lines) and after (solid 
lines) a Tl transformation. 

perpendicularly to its previous orientation, and get a new 
length e' ~ e, which is an arbitrary constant parameter, 
and the connectivity of the four neighboring bubbles is 
modified. This transformation induces a concentration of 
tension forces around this new edge that leads to a rapid 
local relaxation governed by the same dynamical equa- 
tions as the rest of the motion. As the process is very 
local, the 2D description presumably fails and the dissi- 
pation in the two vertical Plateau borders involved in the 
Tl must be of the order of the dissipation on the plates. 
This additional dissipation will be taken into account in 
a future modeling. In any case, if the characteristic time 
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of a relaxation after a Tl remains small in comparison 
with the other time scales, they can be seen as instanta- 
neous processes and their specific dissipation rates do not 
modify the foam behavior. 

4 Large bubble instability : qualitative 
behavior 

4.1 Performed simulations 

We investigated the behavior of a monodisperse foam in 
which a single large bubble has been created. This de- 
fect is produced during the first relaxation step by impos- 
ing a large reference area to an arbitrary bubble chosen 
upstream, called LB. The results presented below were 
obtained with a LB area 20 times larger than the small 
bubbles area. The mean flow velocity is the control pa- 
rameter of our study. Once the foam is equilibrated, the 
flow is turned on. The bubble areas adjust then slightly 
in order to establish the macroscopic pressure gradient 
between upstream and downstream. In the frame of the 
mean flow (in which C is at rest), a small upstream mo- 
tion of the bubbles is then observed, corresponding to an 
area decrease upstream and an area increase downstream. 
Simultaneously, the foam deforms around the large bub- 
ble which tends to move faster than the mean flow. If the 
mean velocity is smaller than a threshold value, a station- 
ary shape is obtained and the flow remains a plug flow. 
Otherwise, Tl's occur around LB which begins to migrate. 

The parameters e, e', hP°, the box size and the time 
step were varied to test the numerical stability. The algo- 
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rithm becomes numerically unstable for a too large time 
step, whose maximal acceptable value is mainly related to 
hP° and to e. The stability domain has not been system- 
atically investigated, but we carefully checked the repro- ^ 
ducibility of our results over a large range of parameters. 
Only e has a noticeable influence on the threshold value 
as physically expected. It may indeed be seen as the char- 
acteristic size of a Plateau border, and is thus related to 
the liquid fraction <f> of the foam by the relation 



L 

40 

x/d 



<P~e 2 d/(hd 2 ) ~ e 2 /(dh), (11) 



Fig. 8. Trajectory of the large bubble center of mass in the 
■mean flow frame in a case of a disordered foam (the flow is in 
with d the typical edge length, e 2 the Plateau border sec- the x direction). The distance between two successive points is 
tion area and h the distance between the plates. Except proportional to its relative velocity. 
for the results shown in Fig. 1111 all results were obtained 

This tip grows by producing few Tl's in the front of the 

with a value of e which corresponds to a liquid fraction of 

large bubble and finally absorbs the full bubble and a new 

0.5% (disregarding the additional water layer wetting the 

cycle begins. This qualitative behavior is in good agree- 

plates). 

ment with experiments P^. This jerky motion is puzzling 
and its extensive study is the aim of a future work. The 

4.2 Disordered and ordered foams 

explanation of the whole instability process requires first 
In a monodisperse disordered foam, we don't observe a to understand the mean behavior of the large bubble. This 
stationary motion of the large bubble for a velocity larger is the aim of the present work. 

than the threshold. The shape of the large bubble, the The following results were thus obtained with an hexag- 
magnitude and the direction of its relative velocity flue- onal network, only distorted by a single defect. A station- 
tuate strongly (Fig [SJ) . When the large bubble moves at ary regime, in which the large bubble moves at constant 
its highest velocity, its shape is elongated in the direction velocity, is reached after a very short transient. Despite 
of its relative motion. The pinning of the large bubble the overall homogeneity of the foam, various kinds of mo- 
by some defect in the disordered foam stops its migra- tions may be observed with various stabilities and proba- 
tion. The bubble shape becomes then elongated in a di- bilities of occurrence. In the most frequent situation, the 
rection perpendicular to that of the mean flow. Finally large bubble migrates solely thanks to Tl's that take place 
this shape destabilizes, by developing a tip downstream, near its front and its rear. The motion is periodic, and the 
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large bubble recovers exactly the same neighborhood after 
a migration of one bubble layer. Only few Tl's, involving 
the first or second LB neighbors, are performed during 
this period. The crystalline organization is perfectly re- 
stored beyond the large bubble. The orientation of the 
relative LB motion may differ strongly from the direction 
of the mean flow. It occurs always along the most favor- 
able crystalline direction (see FigEJ) and the angle between 
the mean flow and the relative velocity of the large bubble 
may thus reach 30°. 

More rarely, for slightly different initial conditions, a 
well defined wake may appear, with few bubbles moving 
behind the large bubble with the same velocity. This is 
observed in experimental flows too. Several choices seem 
thus to be offered to the large bubble, which may explain 
the large fluctuations observed experimentally and numer- 
ically in disordered foams. 

Detailed numerical results found for the velocity thresh- 
old and the relative velocity of the large bubble are pre- 
sented in the next section together with analytical predic- 
tions. They are obtained for LB trajectories in the direc- 
tion of the flow, without wake, which is the most prob- 
able situation. The velocities, averaged over a set of ex- 
periments, are in fair agreement with predicted values as 
shown on Figll4l 

5 Continuous description 

The numerical results are compared hereafter to analyt- 
ical predictions based on a completely different point of 
view. We partially forget the specific foam structure, at 
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Fig. 9. Internal system and external system, separated by the 
dashed circle. 

least outside the large bubble, and we treat it as an elasto- 
plastic continuum medium. This point of view has already 
been successfully adopted by Kabla et al. in a quasi-static 
context |22| . Scaling behaviors are determined to empha- 
size the main physical processes involved, independently 
of any numerical prefactor. We give in the following sub- 
sections the analytical expressions of the various forces 
exerted on the network of small bubbles, per unit surface 
of foam (as seen from above), averaged over a few bubbles. 

5.1 Viscous forces 

Below the threshold, the viscous forces are well deter- 
mined, as the velocity field is v(x, y) = v$\i x everywhere. 
The migration of the large bubble induces of course veloc- 
ity variations in the foam, but the LB may move at high 
speed among the others with only very localized varia- 
tions of vertices (and films) velocities. The migration oc- 
curs through a fracture process : few vertices at the front 
of the large bubble are accelerated until they reach their 
downstream neighbors and undergo a Tl. Then the small 
bubbles initially in contact with these vertices move on the 
sides of LB and their vertices recover the mean velocity 
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vq . The vertices between LB and their new downstream artificial splitting of the total pressure in two terms proved 

neighbors are then accelerated at their turn. A similar to be very powerful, as it allows us to reach a full and sim- 

but inverted transformation occurs behind the large bub- pie analytical description of the complex pressure field, in 

ble. Thus velocities higher than v$ and Tl transformations good agreement with the results of computer simulations, 

are confined in the internal system defined on FigEI In the Hereafter, the equilibrium pressure defined in that way 

external system, the variations of the network deformation will thus be the local reference pressure, remembered to 

due to the LB migration occur without topological rear- be specific to each bubble. Its value is discussed in detail 

rangements and induce only negligible velocity gradients, in paragraph[7J Our method should however be improved 

For these reasons the velocity field is assumed to be to tackle more complex situations in which the system can 

vq\i x in the whole external domain, even above the thresh- only relax through many Tl transformations, 
old. The viscous force per unit foam surface is thus (see The resulting force per unit foam surface due to the 

eqslEJ}, hetero geneous gas pressure is consistently split into two 

p = v{vq)vq u = _ v f y(vo)v x \ ^ terms, namely: 



5.2 Pressure field 



F p = F Pteq - hVP, (13) 



mi C ■ i 4- A- The first term of the right-hand side is by definition ob- 

lne pressure in a toam is heterogeneous, even at equi- ° J 

ti i . .i i.ii 1 1 r , , tained in a foam at equilibrium, as discussed above. It will 
librium. Its value is mainly related to the toam topology 1 ' 

,,,, , r-j c iu 1.1.1 j • play no role in the dynamics. The second term is not very 

and to tne number of sides n ot each bubble and arise 1 J J J 

c ,i i , , . r sensitive to the local structure of the foam, and we assume 

horn the competition between tension forces which tend ' 

i , . i ,i , , , i j c i - i , j . that it is correctly described by the gradient of a smooth 

to shrink the bubbles and pressure forces which tend to J J b 

inflate them. This effect remains present even when the dynamical pressure field defined as P — P P eq . This field 

r a j ■ ,1 • j -,i ij will be derived in section 171 

toam flows and is then superimposed with purely dynam- '—' 

ical pressure variations. To better evidence the dynamical 

contributions, we define the equilibrium part of the pres- 5^3 Tension forces 

sure in a given flowing structure as the pressure in each 

bubble after relaxation to the closest equilibrium shape. Liquid foams have a rather simple behavior at small strain. 

The latter shape is obtained numerically by fixing the line A linear elastic response is obtained in this case with 

£ of vertices and relaxing the rest of the foam. As the foam an clastic coefficient scaling as j/d |251l24| . If the stress 

is only locally deformed, relaxation occurs usually, as de- reaches a given threshold, equally of the order of 7/d, 

sired, with no need of Tl transformation. This somewhat plastic deformation occurs (through Tl's). 
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The distribution of surface forces, due to tension forces, computer simulations, boundary effects on the cell sides 

is first decomposed in the same way as the pressure, lead- are actually disregarded as the problem is solved with peri- 

ing to odic boundary conditions. In contrast, in real experiments 

Ft = F 4)eg + Ft (14) in a Hele-Shaw cell of width 2L there is a frictionless slip 

. , on the lateral smooth boundaries, leading to X v = and 

WltXl J? t,eq = —Ep, e q. 

Considering an elastic and incompressible response of 
the foam we get 



&xy = where a is the elastic stress tensor. Boundary 
conditions have only a logarithmic influence on the force 
fields near the large bubble, varying as ln(D/L). They are 
difficult to determine analytically and are not relevant for 
where X is the displacement vector from the closest equi- our purpose. As an illustration, we nevertheless give here 



F P , eg + ^V 2 X , (15) 



librium shape obtained after relaxing the foam as explained the asymptotic expression of the elastic force exerted on 
previously. The weak foam compressibility is disregarded. thc lar S e bubble with a third condition, X = on the lat- 



eral boundaries, well documented in thc literature in the 
context of viscous hydrodynamics |25U26| . 

The large bubble motion is governed by the resulting 
force induced by the external pressure field and the elas- 
6.1 Driving force tic stress acting on the internal system. This force F x , ori- 



6 Instability threshold and large bubble 
velocity 



ented along the mean flow for symmetry reasons is, noting 
C the boundary of the internal system and n its normal, 



Using the three force expressions eqs. II 211 31151 and ne- 
glecting any inertial term, we get the following equation 

for the external system, F x = h / —Pn x + — 1 ( 2 - x n x + [ - y -\ — ) n y ) dl 



d \ dx \ dx dy 

+ h p] + 2^v 2 x = o . (16) (17) 



The analytical expression of this force has been obtained 

The equation of incomprcssibility div X=0 and thc bound- by analogy with thc pro blem of an infinite cylinder pushed 

ary conditions must be added to eq. M at constant velocity in a viscous fluid between two walls 

The internal system shape is approximated by a disc V which is solved in thc literature 25 . The identification 

of diameter D centered at the LB position. As it tends to of each variable is based on the sim ii arity between thc 

migrate towards positive x, its displacement is 5u x . The equations governing an elastic solid and a viscous fluid, 

parameter 5 remains undetermined at that point, it will As dc t a iled in Appendix \K\ it leads to 
be expressed explicitly in the next section. So we get, as 

first boundary condition for eq. El X = 5u x on T>. In x d d ln(^) — 0.91 + 0.43 + ••• 
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The first term, oriented downstream, is the driving force plied to the internal system, ie that due to the foam films 

for the instability. It is a pressure contribution correspond- around LB, is given by Dt](vo)vq (see eq|5J)- This leads to 

ing to the missing viscous forces in the large bubble. It is the following equation for S, with the use of eq. 1191 
similar to an Archimedes force in which the homogeneous 

gravity field would be replaced by the homogeneous vis- T](vq)voD 2 jhS 
cous force field —r)(vo)vo/du x . The resulting pressure force 



d d 

is proportional to the force field intensity times the area of an< ^' * n ^ e um it D ^> d, 



■V = Dt){v )v , (20) 



the hole. The second contribution, oriented upstream, is ^ t](vq)vqD 2 ^ 
the elastic response of the foam. As expected, it is propor- 
tional to the elastic coefficient 7/c? and to 6. The denomi- The threshold is reached when the maximal stress in 
nator is an infinite sum expanded in the small parameter the foam reaches the yield stress 7/d. The largest elastic 
D/L. Its logarithmic dominating term is specific to the 2D stresses are localized in front of the large bubble and just 
case, in which elastic problems cannot be solved without behind and their values are determined from dimensional 
boundary conditions at finite distance. This explicit ana- arguments. As the elastic coefficient scales as 7/d, they 
lytical solution allows to quantify roughly the influence of scale as jS/(dL geo ), where L geo is a length related to the 
the cell size. In any case, as other numerical prefactors will problem geometry. The box size appears only through log- 
be disregarded in the following, we neglect this logarithmic arithmic corrections. In a purely continuous medium de- 
dependence and we only checked that the denominator re- scription, the small bubble size scales out and thus L geo ~ 
mains close to unity. The driving force expression used in D, which is the value retained in the next part 16.31 Any- 
the next subsection is thus way, the discrete nature of the foam should probably be 
vi v o)voD 2 jhS taken into account here. The small length scale d influ- 

Fx = d T (19) 

ences the radius of curvature of the LB in the region sur- 

From this force, compared to the viscous force exerted on r0U nding the fracture tip, and therefore the stress concen- 

the large bubble itself, we determine the value of 5 in the tration. Both assumptions L geo = D or d are thus investi- 

following. gated below, as each might be valid in some domain which 

depends on the actual aspect ratio. 
6.2 Instability threshold From thiS; we deduce the maximal value S max ~ L geo 



At small velocity, the elastic response of the foam compen- 
sates the driving force and the large bubble has the same 



and the velocity threshold v t h- 



hL 



velocity vq as the rest of the foam. The viscous force ap- ?l( v th)vth ~ 7 ^2 (^) 
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Using T](v t h)vth = l{Vwv/l) a (see eq.EJ we get, with the 
adimensional capillary number Ca = r/ w v/j, 



15 



Ca 



th 



(23) 



V D 2 

Depending on the value of L geo and on the value of 
a (between 0.5 and 0.66 as discussed in section l2~2*|l . the 
model predicts an exponent for D ranging between —1.5 
and —4. The experimental value is of the order of —3.7, 
in a case where a — 0.5 (see Fig llOl from £Q). The choice 
Lgeo — d, which is thus in better agreement with the ex- 
periments than is L geo = D, reinforces the importance of 
the discrete structure of the foam. New numerical simu- 




0.04 0.06 

h/D 



Fig. 10. Experimental results from Cantat et al. '11. Each 
point in the plane (h/D; r/uo/7) represents a flow, with a large 
bubble migrating (circle) or not (square). Two domains ap- 
pear, separated by a boundary correctly fitted by the relation 



lations, especially with a large range of large bubble size, Ca t h = 68(h/D) 3 7 (dashed line). Fluctuations produce 



are nevertheless needed to be fully conclusive. 



lap between the two domains which results in a large uncer- 



The liquid fraction plays an important role in the plas- tainty on the exponent value. The boundary line, obtained from 
tic threshold [231 , and therefore influences the prefactor in eq^with L geo = D and 1/a = 2 and an adjustable multiplica- 
cqE3 In our vertex model simulation, the liquid fraction twe Prefactor, is further plotted (full line). 



is related to the minimal distance e allowed between two 
vertices before a Tl occurs (see ea. Illjl. The variation of 
the threshold with this parameter is displayed on Fig. 1111 
All the simulation results presented on the other graphs 
were performed with a fixed e value, as small as allowed 
by the numerical stability. 

6.3 Large bubble relative velocity 



of the elastic stress. Our model for the LB behavior is thus 
the following : if a Tl occurs at time t — 0, the elastic force 
vanishes, the value of S is zero and the large bubble ve- 
locity is given by eq|20|: ri(v(0))v(0) = rj(vo)voD / d. Then 
S increases with time, inducing an elastic stress increase 
and a LB velocity decrease. Finally i5 reaches its limit D 
at t — t m ax and a new cycle begins (see Fig. [T^l . 



As S(t) = Jq (v — vo)dt the equation of motion becomes, 
For flow velocities larger than the threshold, plastic trans- between t = and t = t max , (see eal2THl 



formations occurs regularly around the large bubble which 
migrates through the foam with a velocity v(t). As these 
Tl's are strongly localized around LB, they can be consid- 



jh/d f\v - v Q )dt + vMvoD 2 = ^ W)u(t) (24) 



As detailed in Appendix [Bl this differential equation 
ered as discrete events relaxing suddenly the largest part is analytically solvable. The approximate solution detailed 
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" numerical results 
— linear fit 



.a 0.9 - 

> 0.8 - 



X 0.5 - 

>° 

> 




E/d 

Fig. 11. Velocity threshold as a function of the minimal dis- 
tance e allowed between two vertices before a Tl occurs. 

below is obtained by neglecting the variation of n(v) with 
v. It differs only by few percents from the exact one and 
is much simpler (see Fig. II 3[). Taking the time derivative, 
we get 



7/1 . . . , dv 



v(t) =v + («(0) - f )e- t/r 

_ dDj] 
7 h 



(25) 

(26) 
(27) 



t/(d/V ) 

Fig. 12. Numerical LB velocity as a function of time, for 
vq > Vth ■ Tl 's involving LB edges occur in avalanches of 2 
or 3 events and are indicated by arrows. At these times the 
LB is suddenly deformed and the center of mass velocity di- 
verges numerically. Just after a Tl avalanche, the elastic stress 
is strongly reduced and the LB velocity is large. As the stress 
accumulates in the small bubbles network the LB velocity de- 
creases until the yield stress is reached and new plastic events 
occur. In contrast, if vo < Vth the LB relative velocity tends to 
zero before the yield stress is reached. 



< V — Vq > 



-Vth/Vo 



(31) 



Expression l26l is valid until <5 reaches 5 max ~ D at time 



\v(t) - v )dt = D = - T (v(0) - v ) (e-*— ' T - l) 

(28) 



So 



t, 



■In [I 



TV 



and finally, with <> denoting the average in time, 

D -D 

< v — vq > 



tmax T l n (l ±- 

TVq 



(29) 



(30) 



voD/d i n ( 1 _sa. 

The velocity threshold is consistently obtained from eq. 
1221 with 77 = const., namely v t h = "fh/r/D. 

This expression is in good agreement with our numer- 
ical simulations (see Fig. I13|l . The small oscillations of 
the numerical points on both sides of the analytical curve 
are related to a bubble structure modification around LB 
at some velocities. The logarithmic behavior above the 
threshold renders the transition almost discontinuous. A 
small modification of the threshold value, due to a small 
polydispersity for example, induces a large variation of 
the LB velocity. The fluctuations observed experimentally 
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A better understanding of fluctuations is beyond the 
scope of this article. It would necessitate an improved de- 
scription of the force field in the foam and of the structural 
disorder. 



x numerical result 
- - analytical result (exact) 
■ ■ analytical result (approx.) - 



mean velocity / threshold velocity 



7 Pressure field. 
7.1 Pressure at equilibrium 



Fig. 13. Large bubble relative velocity as a function of the 

mean flow velocity. The first analytical law is given by eg. \45\ 

and the second one comes from ea. \31\ Numerical vrefactors are m, . . r , , 

1 1 lhe pressure at equilibrium m a 2D loam obeys elegant 



adjusted to obtain the best fit between numerical and analytical 
values. Their values are respectively 0.25 and 0.2, which is of 
the order of unity as expected. The ratio D/d is 3. 



Q 

o 

> 



> 0.02 
> 



□ 0.03 < h/D < 0.04 

♦ 0.04 < h/D < 0.05 

A 0.05 < h/D < 0.07 

t 0.07 < h/D < 0.09 

— analytical result 



_* A 



□ 

-is. A»_ 




0.1 



v /v 

V th 



and simple rules |12|. It is mainly related to the number of 
sides of each bubble. The average number of sides is six as 
a consequence of Euler's relation in 2D The pressure 
inside bubbles with n > 6 is larger than the external ref- 
erence pressure whereas it is smaller when n < 6. When 
bubbles with n ^ 6 can be considered as relatively isolated 
defects in a monodisperse foam, the pressure distribution 
has been calculated by Graner et al. from an electrostatic 
analogy. Bubbles with n ^ 6 play the role of topological 
charges and modify the pressure field in the same way as 
positive or negative charges produce an electric potential. 



Fig. 14. Experimental results from Cantat et al. OF Large Th e induced pressure field around each defect is, with r 



bubble relative velocity as a function of the mean flow velocity, 
rescaled by the threshold velocity obtained on Fia MtA for var- 
ious large bubble sizes. The analytical fit comes from eq. 1311 
with a single adjustable prefactor. 



the distance to the defect, (see |27p 



P P° = -l^tlnC-) 
d 3 V 



(32) 



The logarithmic divergence is specific to the 2D geome- 
try and we recover an expression similar to the electric 
are therefore relatively important, but the mean veloci- potential produced by an infinite and uniformly charged 
ties, averaged over many experiments, are very well fitted line. The various contributions are additive, and a mul- 
by the previous law as shown on Fig^J from Cantat and tipolar expansion can thus be used to derive the pres- 



Delannay 0. 



sure field far from the defects. During its migration, the 
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large bubble adopts an elongated shape which is kept even 
after relaxation. The distribution of 5-sided bubbles on 
both sides of the large bubble with n^B sides is equiva- 
lent to a quadrupolar distribution of topological charges, 
from which we deduce the dominant term in the pressure 
variation, up to a numerical prefactor, 

2y 2 x 2 



(P - P°) ~ n LB d"f 



(33) 



(x 2 + y 2 ) 2 (x 2 + y 2 ) 2 j 

The results of numerical simulations are in good agree- 
ment with the predictions of eq. 1331 as shown in figure 
Figdi-b. 

7.2 Pressure during the flow 

The pressure field in the flowing foam has already been 
discussed in detail in 8 , and we just recall the main re- 
sults here. As it is less sensitive to long range effects than 
the displacement field, it can be computed analytically 
without taking into account the positions of the bound- 
aries. The final equation for the pressure is , with ro the 
position of the large bubble, 

TJVqX T/VqD 2 X — Xq 



p = — 



(34) 



d 2nd (r-r ) 2 ' 

The first term is the linear pressure variation responsi- 
ble for the flow. The second term, compared to numerical 
results on Fig llfil is related to the foam deformation. 

In contrast to the deformation of the foam and with 
the resulting tension forces, the pressure cannot be mea- 
sured locally from images of a 2D foam. As shown on Fig. 
1161 it may nevertheless be very heterogeneous and plays 
a crucial role in the foam dynamics. It is thus an experi- 





• numerical results 
— y=-4/x 2 



30 35 



Fig. 15. Pressure field around the large bubble, obtained by 
relaxation from an out of equilibrium situation, with a velocity 
larger than the threshold value, (a) Pressure field in gray level : 
pressure decreases from white to black (for clarity, the color of 
LB is arbitrarily chosen, its real pressure being lower than that 
of small bubbles), (b) Pressure along the horizontal segment 
y — from LB to the right boundary. The solid line is a fit 
to the analytical expression given by ea \3'J\ with an adjustable 
prefactor. 

mental challenge to determine the pressure in each bubble 
in a flowing foam. 

The same approach can in principle be used to measure 
the displacement field. We thus compared the position of 
the bubble center before and after relaxation. The small 
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ment we want to study. If this contribution is subtracted, 
the displacement field is nevertheless qualitatively differ- 
ent from the purely elastic displacement field computed 
with the software Freefem, for the boundary conditions 
obtained with the vertex model simulations. More simula- 
tions are needed to conclude to the physical relevance of 
such results and to exclude that they are just artefacts. 
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Fig. 16. Pressure field in the foam after removing the equi- 
librium and linear contributions. The foam flows from the left 
(x < 0) to the right {x > 0). (a) Numerical results. Light gray 
corresponds to high pressure values, (b) Another view of the 
same results : pressure in the foam along the line y — const., 
for three value of y. The line y=0 goes from upstream to down- 
stream, trough the middle of LB. The values are well fitted by 
the analytical expression \y4\ 

compressibility of the foam induces a mean parabolic dis- 
placement of the form X x = A(x 2 — L 2 ) because the bubble 
volumes increase upstream and decrease downstream dur- 
ing relaxation. Even with very small variations of bubble 
areas, as the displacement is the integral of the deforma- 
tion, this effect is larger than the small elastic displace- 



8 Summary and conclusion 

In conclusion, we have shown that the large bubble insta- 
bility can be reproduced with simple numerical simula- 
tions, keeping only the skeleton of the foam structure and 
the three key ingredients which are tension forces, vis- 
cous forces and volume conservation in the bubbles. The 
main physical aspects of this dynamical behavior is under- 
stood as a complex interplay between the viscous, elastic 
and plastic properties of the foam, seen as a continuum 
medium. This theoretical approach leads to a good agree- 
ment with numerical and experimental results. It can be 
adapted for various situations, as foam flows around ob- 
stacles for instance. Future work will have to take into 
account the disordered (or possibly crystallized) foam or- 
ganization to deal with the velocity fluctuations of the 
large bubble. 

A very important perspective is to predict the flow 
properties of a fully polydisperse foam. A mean field the- 
ory will probably be insufficient, as bubbles seem to be 
strongly coupled and as the flowing foam reorganizes it- 
self with a resulting size segregation effect. 
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A Hydrodynamical analogy 

Hydrodynamics equations for an incompressible viscous 
fluid of velocity Vh, pressure P^ and viscosity 77 is : 

r,Av h - VP h = ; V • v h = . (35) 



Eauation l24l is adimensional with the time unit r = D d Ca'J 
With, t = t/r and v = v/vq and omitting the bar notation, 
we get, 



t(v - 1)dt = 1 - r 2/3 (39) 



whose differentiation yields 



^4/3 _ v l/3 + f| = (40) 



The solution of this differential equation with separable 
The problem of an infinite cylinder of diameter D pushed variables is 

at constant velocity Vh between two walls separated by t = H(v(t)) — H(v(0)) (41) 

a distance 2L has been solved by Faxen [21]. The vari- 



with 



ables can be identified as follows 77V/, <-► jh/d X, Pj 



1, 



( VxW + xW + l \ ^ _ x ( 2x^ + 1 
rjvox/d + P and Vk <-> S. The expression of the force ex- n W = m \ x i/ 3 _ 1 I _v<3 tan ^ 

erted on the cylinder is |25U2(ij 

(42) 

The initial velocity is v (0) = (D/d fl 2 and the maximum 

-Ph,n x + T) \2^^n x + (^IH. + \ j s 



q x x \dx dyJ J time before a Tl is given by 

' ln{ lL) _ 0.9157 + 1.73 ^ + ... ' " ' ' f m {v _ l)dt |«(t m ) 2 /3 = 2_ (43) 

Jo d 2 v t 



2 ( D D x 3/2 



Using the variable identification, we obtain 

-(r,v x/d + P)n x + " (tm)= 3U"^j ^ 

I dX y dX x \ \ fj na Uy express the mean relative velocity of the large 



7 / , f ox x fdx y dx x \ 



Air^hS bubble in terms of adimensional variables as: 

(37) 



d (ln(%)- 0.9 + 1.73^ + 

The term we need to compute is given by eq. 1171 and is 
thus 



F^/^fdx-—^ ^ n2 (38) 



ln{%- 0.9 + 1-73^ + 



D 1 

< v ~ 1 >_ ^ffwa - HWJdW) 
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